      Program equilib
      !****************************************************************
      include 'PPMA.Fort'
      include 'STATD.Fort'
      Dimension r(0:IMax),alpha(0:IMax),
     +		theta(0:JMax),
     +          phi(0:IMax,0:JMax),
     +          rho(0:I0,0:JMax)
      real etime,dtime,tarray(2)
      Common /pi/pi
      Logical Log
      Logical LogInit
      Logical LogSt
      data LogSt/.true./
      Dimension r2(0:10)

      i=errtrap('4'x) ! Convex: trap at math. errors
      r2(0)=1d0
      r2(1)=0.9d0
      r2(2)=0.8d0
      r2(3)=0.7d0
      r2(4)=0.65d0
      r2(5)=0.625d0
      r2(6)=0.620d0
      r2(7)=0.615d0
      r2(8)=0.610d0
      r2(9)=0.600d0
      open(1,file='output')
      CN=1.5d0
      pi=dacos(-one)
      write(*,*) pi,dcos(pi)
      dr=one/I0
      do i=0,I0
          r(i)=dr*i
      enddo
      dalpha=one/r(I0-1)-one/r(I0)
      do l=0,IMax-I0+1
          alpha(l)=dalpha*l
      enddo
      do i=I0+1,IMax-1
          l=IMax-i
          r(i)=one/alpha(l)
      enddo
      dtheta=half*pi/JMax
      do j=0,JMax
	  theta(j)=dtheta*j
      enddo
      i1=i0
      j1=JMax
c     === initial phi ==>
      if(LogSt) then ! === spheroidal structure ==>
	  i20=i0
          j2=0
          do j=0,JMax
              do i=0,IMax
                  if(i.le.i20) then
                      phi(i,j)=-1.5d0
                      if(i.eq.i0) phi(i,j)=-0.5d0
                  else
                      phi(i,j)=-0.1d0
                  endif
              enddo
          enddo
      else ! === ring-like structure ==>
	  i20=i0/2
	  j2=JMax
	  Bs=-one
	  cpsi0=0.2d0
	  phi1=Bs-cpsi0*Fpsi0(r(i1),theta(j1))
	  phi2=Bs-cpsi0*Fpsi0(r(i20),theta(j2))
          do j=0,JMax
              do i=0,IMax
                  if(i.le.i20) then
		      phi(i,j)=phi2
                  else
                      if(i.le.i1) then
			  if(i.ne.i1) then
	                      phi(i,j)=10d0*(Bs-cpsi0*Fpsi0(r(i),theta(j2)))
			  else
                              phi(i,j)=phi1
			  endif
	              else
                          phi(i,j)=-0.1d0
		      endif
		  endif
              enddo
          enddo
      endif
c     <== initial phi ===
      write(1,*) 'time',etime(tarray),dtime(tarray)
      do k=0,9
          if(LogSt) then ! === spheroidal structure ==>
              i2=r2(k)*i0
	      iH=0
	      jH=JMax
          else ! === ring-like structure ==>
	      i2=i20+k
	      iH=(i1+i2)/2
	      jH=JMax
          endif
	  LogInit=.true.
	  write(*,*) 'CN,LogSt,i2', CN,LogSt,i2
	  write(1,*) 'CN,LogSt,i2', CN,LogSt,i2
          Bs=( phi(i2,j2)*Fpsi0(r(i1),theta(j1))
     +        -phi(i1,j1)*Fpsi0(r(i2),theta(j2)))
     +       /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
          idelta=I0/2
          write(1,*) 'n,r'
          write(1,1) n,(r(i),i=0,IMax,idelta)
          write(1,*) 'n,phi,Bs'
	  n_max=60
	  log=.true.
	  n=1
          do while(log.and.(n.le.n_max))
	      Bs0=Bs
              Call new_phi(phi,r,dr,alpha,dalpha,theta,dtheta,Bs,CH,CN,
     +                     i1,j1,i2,j2,iH,jH,LogInit)
	      if(dabs((Bs0-Bs)/Bs).lt.1d-7) log=.false.
              if(mod(n, 1).eq.0) write(1,1) n,(phi(i,0),i=0,IMax,idelta),Bs
              write(*,*) n,phi(iH,jH),Bs
	      n=n+1
          enddo
1         Format(i5,4(1pe10.2),1pe28.20)
          write(1,*) 'time',etime(tarray),dtime(tarray)
          cpsi0=(phi(i2,j2)-phi(i1,j1))
     +          /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
c          write(1,*) 'j, phi'
c	  do j=0,JMax
c              write(1,1) j,(phi(i,j),i=0,IMax,idelta)
c          enddo
          Call Integrals(CM,V,CJ,T,W,Pip,phi,r,theta,dtheta,
     +                   i1,j1,i2,j2,iH,jH,Bs,CH,CN)
	  VT=abs(2d0*T+W+3d0*Pip)/abs(W)
          pmax=zero
          do i=0,I0
              do j=0,JMax
                  psi=cpsi0*Fpsi0(r(i),theta(j))
                  rho(i,j)=FRho(psi,phi(i,j),Bs,CH,CN)
		  pmax=max(pmax,Fp(rho(i,j),psi,phi(i,j),Bs,CH,CN))
              enddo
          enddo
	  write(1,*) '  k  ',
     +               '   r(i2)   ',
     +               '   cpsi0   ',
     +               '    CM     ',
     +               '     V     ',
     +               '    CJ     '
	  write(1,2) k,r(i2),cpsi0,CM,V,CJ
	  write(1,*) '     ',
     +               '     T     ',
     +               '     W     ',
     +               '    VT     ',
     +               '   pmax    ',
     +               '  T/abs(W) '
	  write(1,3) T,W,VT,pmax,T/abs(W)
2         Format(i5,6(1pe11.3))
3         Format(5x,6(1pe11.3))
c          write(1,*) 'rho'
c          Call Plot(rho,0,I0,0,JMax,1)
      enddo
      write(1,*) 'time',etime(tarray),dtime(tarray)
      stop
      end

      Subroutine new_phi(phi,r,dr,alpha,dalpha,theta,dtheta,Bs,CH,
     + 			 CN,i1,j1,i2,j2,iH,jH,LogInit)
      !****************************************************************
      include 'PPMA.Fort'
      include 'STATD.Fort'
      Parameter(KMax=(IMax+1)*(JMax+1))
      Parameter(nnzero=8*KMax)
c      Parameter(LICN=3*nnzero,LIRN=3*nnzero/2)
      Parameter(LICN=20*nnzero,LIRN=20*nnzero/2)
      Parameter(NIN=5,NOUT=1)
      Dimension r(0:IMax),alpha(0:IMax),theta(0:JMax),
     +          phi(0:IMax,0:JMax),
     +          rho(0:IMax,0:JMax),
     +          a_1(0:IMax),b_1(0:IMax),c_1(0:IMax),
     +          a_2(0:JMax),b_2(0:JMax),c_2(0:JMax),
     +          b(KMax),
     +          bim(KMax),
     +          bip(KMax),
     +          bjm(KMax),
     +          bjp(KMax),
     +          bkH(KMax),
     +          bk1(KMax),
     +          bk2(KMax),
     +          b0(0:JMax),
     +          ff(KMax)
      Dimension ir0(nnzero),ic0(nnzero)
      Dimension A(LICN),
     +          W(KMax),
     +          ICN(LICN),
     +          IDISP(10),
     +          IKEEP(KMax,5),
     +          IRN(LIRN),
     +          IW(KMax,8),
     +          IVECT(nnzero),
     +          JVECT(nnzero)
      Logical GROW,LBLOCK,ABORT(4)
      Logical LogInit
      Common /pi/pi

      k1=kk(i1,j1)
      k2=kk(i2,j2)
      kH=kk(iH,jH)
      cpsi0=(phi(i2,j2)-phi(i1,j1))
     +      /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
      psi_iH_jH=cpsi0*Fpsi0(r(iH),theta(jH))
      CH=Bs-phi(iH,jH)-psi_iH_jH
      do i=1,KMax
          b(i)=zero
          bim(i)=zero
          bip(i)=zero
          bjm(i)=zero
          bjp(i)=zero
          bkH(i)=zero
          bk1(i)=zero
          bk2(i)=zero
          ff(i)=zero
      enddo
      do i=1,I0-1
          a_1(i)=one/(r(i)*dr)**2*(half*(r(i-1)+r(i)))**2
          c_1(i)=one/(r(i)*dr)**2*(half*(r(i)+r(i+1)))**2
          b_1(i)=-a_1(i)-c_1(i)
      enddo
      do l=1,IMax-I0
          i=IMax-l
          a_1(i)=alpha(l)**4/dalpha**2
          c_1(i)=alpha(l)**4/dalpha**2
          b_1(i)=-a_1(i)-c_1(i)
      enddo
      do j=0,JMax
          if((j.eq.0).or.(j.eq.JMax)) then
              if(j.eq.0) then
                  c_2(j)=two**2/dtheta**2
                  b_2(j)=-c_2(j)
              else
                  a_2(j)=two/dtheta**2
                  b_2(j)=-a_2(j)
              endif
          else
              a_2(j)=one/(sin(theta(j))*dtheta**2)
     +               *half*(sin(theta(j-1))+sin(theta(j)))
              c_2(j)=one/(sin(theta(j))*dtheta**2)
     +               *half*(sin(theta(j))+sin(theta(j+1)))
              b_2(j)=-a_2(j)-c_2(j)
          endif
      enddo
c     === aproximation of delta-operator ==>
c     === internal space ==>
      do i=1,IMax-1
          do j=0,JMax
	      k  =kk(i,j)
	      b(k)=b_1(i)+b_2(j)/r(i)**2
	      bim(k)=a_1(i)
              bip(k)=c_1(i)
	      if(j.ne.0) then
                  bjm(k)=a_2(j)/r(i)**2
              endif
	      if(j.ne.JMax) then
		  bjp(k)=c_2(j)/r(i)**2
              endif
          enddo
      enddo
c     <== internal space ===
c     === r=0 ==>
      b(kk(0,JMax))=-6d0/dr**2
      do j=0,JMax
          theta1=half
     +           *(theta(max(0,j-1))+theta(j))
          theta2=half
     +           *(theta(j)+theta(min(JMax,j+1)))
	  b0(j)=6d0/dr**2*(cos(theta1)-cos(theta2))
      enddo
      do j=0,JMax-1
	  k=kk(0,j)
	  b(k)=one
	  bjp(k)=-one
      enddo
c     <== r=0 ===
c     === alpha=0 ==>
      do j=0,JMax
	  b(kk(IMax,j))=one
      enddo
c     <== alpha=0 ===
c     <== aproximation of delta-operator ===
      do i=0,IMax-1
          do j=0,JMax
	    k=kk(i,j)
	    if(k.ge.kk(0,JMax)) then
	      psi=cpsi0*Fpsi0(r(i),theta(j))
              ff(k)=4d0*pi
     +              *( FRho(psi,phi(i,j),Bs,CH,CN)
     +                +FdRho_dphi(psi,phi(i,j),Bs,CH,CN)
     +                 *(-phi(i,j)+Bs-psi)
     +                -FdRho_dCH(psi,phi(i,j),Bs,CH,CN)*CH
     +               )
	      b(k)=b(k)-4d0*pi*FdRho_dphi(psi,phi(i,j),Bs,CH,CN)
              bk1(k)= bk1(k)
     +               -4d0*pi*FdRho_dphi(psi,phi(i,j),Bs,CH,CN)
     +                *Fpsi0(r(i2),theta(j2))
     +                /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
              bk2(k)= bk2(k)
     +               +4d0*pi*FdRho_dphi(psi,phi(i,j),Bs,CH,CN)
     +                *Fpsi0(r(i1),theta(j1))
     +                /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
              bk1(k)= bk1(k)
     +               +4d0*pi*FdRho_dphi(psi,phi(i,j),Bs,CH,CN)
     +                *Fpsi0(r(i),theta(j))
     +                /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
              bk2(k)= bk2(k)
     +               -4d0*pi*FdRho_dphi(psi,phi(i,j),Bs,CH,CN)
     +                *Fpsi0(r(i),theta(j))
     +                /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
              bkH(k)= bkH(k)
     +               +4d0*pi*FdRho_dCH(psi,phi(i,j),Bs,CH,CN)
              bkH(k)= bkH(k)
     +               +4d0*pi*FdRho_dCH(psi,phi(i,j),Bs,CH,CN)
     +                *Fpsi0(r(i2),theta(j2))
     +                /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
              bk2(k)= bk2(k)
     +               -4d0*pi*FdRho_dCH(psi,phi(i,j),Bs,CH,CN)
     +                *Fpsi0(r(i1),theta(j1))
     +                /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
              bk1(k)= bk1(k)
     +               -4d0*pi*FdRho_dCH(psi,phi(i,j),Bs,CH,CN)
     +                *Fpsi0(r(iH),theta(jH))
     +                  /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
              bk2(k)= bk2(k)
     +               +4d0*pi*FdRho_dCH(psi,phi(i,j),Bs,CH,CN)
     +                *Fpsi0(r(iH),theta(jH))
     +                /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
	    endif
          enddo
      enddo
c     === forming initial data for NAG program ==>
      NZ=0
      do k=1,KMax
          i=(k-1)/(JMax+1)
	  j=k-1-i*(JMax+1)
          NZ=NZ+1
	  A(NZ)=b(k) ! main diagonal
	  ir0(NZ)=k
	  ic0(NZ)=k
          if(i.ne.0) then
              NZ=NZ+1
	      A(NZ)=bim(k) ! papall. to main diagonal
	      ir0(NZ)=k
	      ic0(NZ)=kk(i-1,j)
          endif
          if(i.ne.IMax) then
              NZ=NZ+1
	      A(NZ)=bip(k) ! papall. to main diagonal
	      ir0(NZ)=k
	      ic0(NZ)=kk(i+1,j)
          endif
          if(j.ne.0) then
              NZ=NZ+1
	      A(NZ)=bjm(k) ! papall. to main diagonal
	      ir0(NZ)=k
	      ic0(NZ)=kk(i,j-1)
          endif
          if(j.ne.JMax) then
              NZ=NZ+1
	      A(NZ)=bjp(k) ! papall. to main diagonal
	      ir0(NZ)=k
	      ic0(NZ)=kk(i,j+1)
          endif
          NZ=NZ+1
	  A(NZ)=bkH(k) ! point 0
	  ir0(NZ)=k
	  ic0(NZ)=kH
          NZ=NZ+1
	  A(NZ)=bk1(k) ! point 1
	  ir0(NZ)=k
	  ic0(NZ)=k1
          NZ=NZ+1
	  A(NZ)=bk2(k) ! point 2
	  ir0(NZ)=k
	  ic0(NZ)=k2
      enddo
      do j=0,JMax
          NZ=NZ+1
          A(NZ)=b0(j) ! \Delta\Phi at r->0
	  ir0(NZ)=kk(0,JMax)
	  ic0(NZ)=kk(1,j)
      enddo
      do n=1,NZ
          if(LogInit) then
	      IRN(n)=ir0(n)
	      ICN(n)=ic0(n)
          else
	      IVECT(n)=ir0(n)
              JVECT(n)=ic0(n)
          endif
      enddo
      U=0.1d0
      LBLOCK=.TRUE.
      GROW=.TRUE.
      ABORT(1)=.TRUE.
      ABORT(2)=.TRUE.
      ABORT(3)=.FALSE.
      ABORT(4)=.FALSE.
      ! IFAIL=110
      IFAIL=010
      CALL X04ABF(1,NOUT)
      if(LogInit) then
c          CALL F01BRF(KMax,NZ,A,LICN,IRN,LIRN,ICN,U,IKEEP,IW,W,LBLOCK,GROW,
c     + 	              ABORT,IDISP,IFAIL)
c	  THR=1d-2
	  THR=1d-8
          CALL M28Y12(KMax,NZ,A,LICN,IRN,LIRN,ICN,U,IKEEP,IW,W,LBLOCK,GROW,
     + 	              ABORT,IDISP,IFAIL,THR)
          LogInit=.false.
          LogInit=.true. ! only M28BYS (25 October)
	  NM28BYS=0
      else
          ETA=1d-4
          ! IFAIL=110
          IFAIL=010
c          CALL F01BSF(KMax,NZ,A,LICN,IVECT,JVECT,ICN,IKEEP,IW,W,GROW,ETA,
c     +                RPMIN,ABORT(4),IDISP,IFAIL)
          CALL M28BYS(KMax,NZ,A,LICN,IVECT,JVECT,ICN,IKEEP,IW,W,GROW,ETA,
     +                RPMIN,ABORT(4),IDISP,IFAIL)
	  NM28BYS=NM28BYS+1
	  if(NM28BYS.ge.2) LogInit=.true.
      endif
      MTYPE=1
      CALL F04AXF(KMax,A,LICN,ICN,IKEEP,ff,W,MTYPE,IDISP,RESID)
c     <== forming of initial data for NAG program ===
      do k=1,KMax
          i=(k-1)/(JMax+1)
	  j=k-1-i*(JMax+1)
          phi(i,j)=ff(k)
      enddo
      Bs=( phi(i2,j2)*Fpsi0(r(i1),theta(j1))
     +    -phi(i1,j1)*Fpsi0(r(i2),theta(j2)))
     +   /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
      cpsi0=(phi(i2,j2)-phi(i1,j1))
     +      /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
      psi_iH_jH=cpsi0*Fpsi0(r(iH),theta(jH))
      CH=Bs-phi(iH,jH)-psi_iH_jH
      return
      end

      integer Function kk(i,j)
      !****************************************************************
      include 'STATD.Fort'
      kk=1+i*(JMax+1)+j
      return
      end

      double precision Function FRho(psi,phi,Bs,CH,CN)
      !****************************************************************
      include 'PPMA.Fort'
      if(CN.eq.zero) then
          if((Bs-phi-psi).gt.zero) then
              FRho=one
          else
              FRho=zero
          endif
      else
          if((Bs-phi-psi).gt.zero) then
              FRho=((Bs-phi-psi)/CH)**CN
          else
              FRho=zero
          endif
      endif
      return
      end

      double precision Function Fp(rho,psi,phi,Bs,CH,CN)
      !****************************************************************
      include 'PPMA.Fort'
      if(CN.eq.zero) then
          Fp=Bs-phi-psi
      else
	  Fp=(CH/(CN+one))*rho**(one+one/CN)
      endif
      return
      end

      double precision Function FdRho_dphi(psi,phi,Bs,CH,CN)
      !****************************************************************
      include 'PPMA.Fort'
      if(CN.eq.zero) then
          FdRho_dphi=zero
      else
          if((Bs-phi-psi).gt.zero) then
              FdRho_dphi=-CN*(Bs-phi-psi)**(CN-1)/CH**CN
          else
              FdRho_dphi=zero
          endif
      endif
      return
      end

      double precision Function FdRho_dCH(psi,phi,Bs,CH,CN)
      !****************************************************************
      include 'PPMA.Fort'
      if(CN.eq.zero) then
          FdRho_dCH=zero
      else
          if((Bs-phi-psi).gt.zero) then
              FdRho_dCH=(Bs-phi-psi)**CN*(-CN)/CH**(CN+1)
          else
              FdRho_dCH=zero
          endif
      endif
      return
      end

      double precision Function Fpsi0(r,theta)
      !****************************************************************
      include 'PPMA.Fort'
      if(r*sin(theta).le.one) then
          Fpsi0=-half*r**2*sin(theta)**2
      else
          Fpsi0=-half
      endif
      return
      end

      double precision Function Fdpsi_domega(omega)
      !****************************************************************
      include 'PPMA.Fort'
      if(omega.le.one) then
          Fdpsi_domega=-omega
      else
          Fdpsi_domega=zero
      endif
      return
      end

      Subroutine Integrals(CM,V,CJ,T,W,Pip,phi,r,theta,dtheta,
     +                     i1,j1,i2,j2,iH,jH,Bs,CH,CN)
      !****************************************************************
      include 'PPMA.Fort'
      include 'STATD.Fort'
      Dimension r(0:IMax),theta(0:JMax),
     +          phi(0:IMax,0:JMax)
      Common /pi/pi

      cpsi0=(phi(i2,j2)-phi(i1,j1))
     +      /(Fpsi0(r(i1),theta(j1))-Fpsi0(r(i2),theta(j2)))
      psi_iH_jH=cpsi0*Fpsi0(r(iH),theta(jH))
c      write(*,*) 'cpsi0',cpsi0
      CH=Bs-phi(iH,jH)-psi_iH_jH
      CM=zero
      V=zero
      CJ=zero
      T=zero
      W=zero
      Pip=zero
      do i=0,i0
          do j=0,JMax
              psi=cpsi0*Fpsi0(r(i),theta(j))
              rho=FRho(psi,phi(i,j),Bs,CH,CN)
     	      c1=-cpsi0
     +            *Fdpsi_domega(r(i)*sin(theta(j)))
     +            *r(i)*sin(theta(j))
	      vphi=sqrt(max(zero,c1))
	      rm=half*(r(max(0,i-1))+r(i))
	      rp=half*(r(i)+r(i+1))
	      thetam=half*(theta(max(0,j-1))+theta(j))
	      thetap=half*(theta(j)+theta(min(JMax,j+1)))
              dV=(rp**3-rm**3)/3d0
     +           *(cos(thetam)-cos(thetap))
     +           *4d0*pi
	      CM=CM+rho*dV
	      if(rho.ge.1d-15) then
		  V=V+dV
	          Pip=Pip+Fp(rho,psi,phi(i,j),Bs,CH,CN)*dV
	      endif
	      CJ= CJ
     +		 +rho*vphi
     +            *r(i)*sin(theta(j))
     +            *dV
	      T=T+rho*vphi**2/2d0*dV
	      W=W+rho*phi(i,j)/2d0*dV
          enddo
      enddo
      return
      end

      Subroutine Plot(UU,l1,l2,m1,m2,NCh)
      !****************************************************************
      ! plot UU(*,*)
      include 'PPMA.Fort'
      Dimension UU(L1:L2,M1:M2)
      Character G(0:35)

      G( 0) = ' '
      G( 1) = '1'
      G( 2) = '2'
      G( 3) = '3'
      G( 4) = '4'
      G( 5) = '5'
      G( 6) = '6'
      G( 7) = '7'
      G( 8) = '8'
      G( 9) = '9'
      G(10) = 'A'
      G(11) = 'B'
      G(12) = 'C'
      G(13) = 'D'
      G(14) = 'E'
      G(15) = 'F'
      G(16) = 'G'
      G(17) = 'H'
      G(18) = 'I'
      G(19) = 'J'
      G(20) = 'K'
      G(21) = 'L'
      G(22) = 'M'
      G(23) = 'N'
      G(24) = 'O'
      G(25) = 'P'
      G(26) = 'Q'
      G(27) = 'R'
      G(28) = 'S'
      G(29) = 'T'
      G(30) = 'U'
      G(31) = 'V'
      G(32) = 'W'
      G(33) = 'X'
      G(34) = 'Y'
      G(35) = 'Z'
      UUMin=UU(L1,M1)
      UUMax=UU(L1,M1)
      do l=L1,L2
          do m=M1,M2
              UUMin=Min(UUMin,UU(l,m))
              UUMax=Max(UUMax,UU(l,m))
          enddo
      enddo
      Write(NCh,*) ' Plot :  UUMin,UUMax =',UUMin,UUMax
      if(UUMin.lt.UUMax) then
*         do i1=0,8
*             write(NCh,4) (G(i1+i2*9),UUMin+(UUMax-UUMin)/float(35)
*    +                                       *(i1+i2*9),i2=0,3)
*         enddo
*4         Format(4('   ',A1,'->',1pe9.2))
          do m=M2,M1,-1
              write(NCh,1) (G((UU(l,m)-UUMin)
     +                      /(UUMax-UUMin)*35+one/two),l=L1,L2)
          enddo
1         Format(150A1)
      endif
      return
      end
